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INTRODUCTION 

Bubbles are essential particles in many industrial as well 
as natural processes. Heat transfer through boiling, fluid cavi- 
tation, bubble driven circulation systems in metal processing, 
and bubble/free surface interactions in oceans are just a few 
examples of the many roles that bubbles play in the physical 
systems. One of the key questions in phase change dynamics 
is the prediction of mass transfer rate and its effects on the 
fluid flow and heat transfer of the system. There has been 
extensive investigations about bubble growth rate in the last 
few decades. The early work was concerned about solving 
separately the momentum and heat conduction equation and 
can be divided into two main categories: growth rate con- 
trolled by inertia and growth rate controlled by heat diffusion. 
In the inertially-controlled growth the temperature difference 
between bubble and liquid is negligible and the growth is con- 
trolled by the pressure difference inside and outside the bubble 
and inertia of the liquid. Rayleigh (1917) made such assump- 
tion to solve for the collapse of an empty cavity in a large mass 
of liquid. Rayleigh’s solution was confined to spherically 
symmetric, isothermal, and inviscid flow. He also neglected 
the surface tension. Later on, his analysis was extended to 
include liquid viscosity and surface tension. In the heat dif- 
fusion controlled growth the pressure difference inside and 
outside bubble is negligible and the growth is controlled by 
heat transfer from the liquid to the bubble. The investigations 
of Bosnjakovic (1930), Plesset & Zwick (1952), and Skinner 
& Bankoff (1965) are a few examples of analyses based on the 
heat diffusion controlled growth. Subsequent analyses to solve 
the coupled equations have produced several asymptotic solu- 
tions applicable to long time only and a number of complete 
but approximate solutions involving simplifying assumptions 
about the heat conduction equation. See, for example, Plesset 
& Zwick (1954), Scrivcn (1959), and Bankoff (1958) for the 
former and Kosky (1968), Mikic et al. (1970), and Theofanous 
& Patel (1975) for the latter. 

The use of numerical modeling in solving fluid flow with 
phase change is relatively new. One of the earliest of such 
studies is the potential flow solution of Plesset & Chapman 
(1971) for the collapse of an isothermal spherical vapor bub- 
ble in the neighborhood of a solid wall. A number of authors 
have tried to solve the boiling bubble problem by coupling 
the momentum equation with a simplified form of the energy 
equation. See Theofanous et al. (1969) and Prosperetti & 
Plesset (1977), for example. Dalle Donne & Ferranti (1975) 


were among the first to solve the complete equations of motion 
and energy. More recent numerical studies about boiling bub- 
bles include boundary fitted/finite element method of Schunk 
& Rao (1994), finite volume/moving mesh method of Welch 
(1995), level set method of Son & Dhir (1 997), and front track- 
ing/finite difference technique of Juric & Tryggvason (1995). 

So far, the last two studies present the most comprehensive 
numerical modeling of the phase change problem leaving aside 
their two-dimensionality. While it is possible to predict some 
of the qualitative features of the phase change phenomena 
by two-dimensional simulations, it is generally believed that 
these computations fail to provide an accurate account of the 
quantitative features. Therefore, it is necessary to relax this 
restriction in order to obtain a more realistic picture of the 
problem. Our aim is to present a numerical technique for 
modeling three-dimensional fluid flow with phase change. To 
the best of our knowledge, this is the first numerical study 
that addresses this issue. In an earlier study (i.e. Esmaeeli 
& Arpaci, 1998), we used a similar methodology to simulate 
phase change without fluid flow. 

FORMULATION AND NUMERICAL METHOD 

Consider a domain consisting of a liquid and its vapor un- 
dergoing phase change. The material properties of the fluids 
are different but constant within each phase. The governing 
equations for this problem arc conservation of mass, momen- 
tum, and energy equations within each phase, and the jump 
conditions across the interface. Rather than writing the gov- 
erning equations separately for each of the fluids, a single set 
of equations arc used which arc valid for the entire flow field 
and take the jump in properties across the interface into ac- 
count. The momentum equation in conservative form for such 
a flow is: 

^ + V p«¥ = -V/) + V /.(V« + V17 t ) (1) 

dt 

—P9 + j> <THnS{x — T jr )dA. 

In the above equation « is the velocity, p is the pressure, p 
is the density, p is the viscosity, 7j is the gravity, <r is the 
surface tension coefficient, k is twice the mean curvature, and 
W is the unit vector normal to the interface. £(j - T f ) is a 
delta function which is zxro everywhere except at the interface 
where x = i f and dA is differential area of the interface. 
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The energy equation in conservative form is: 

+ v ■ pcirr = v • *vr + (2) 

dt 

<j) hpf(lii — uj).nS(x — X* )dA t 

where T is the temperature, r is the heat capacity, and k is the 
heat conductivity. The last term in the above equation acts as 
a heat source (sink) which is zero away from the interface and 
its inclusion enforces the conventional energy jump condition 
across the front. Here, pf and ¥/ are the fluid density and 
velocity at the liquid or vapor side of the interface, and u ? 
is the interface velocity, h is a measure of difference in the 
enthalpies of the liquid and the vapor and is derived using 
thermodynamics consideration: 

h = h jg + Tie q( ci — c v ). (3) 


Here, hj y is the latent heat of evaporation, T eq is the equilib- 
rium temperature, and ci and r t are the heat capacities of the 
liquid and the vapor, respectively. 

The mass conservation is: 

V • pit = (j) (p v — pi)u>n&(* ~~ 4, W 

where, u n — The above equation is equivalent to the 

conventional mass conservation equation but is better suited 
for numerical formulation. 

The introduction of the interface velocity u t adds an ad- 
ditional unknown to the problem. This new unknown is be- 
ing taken care of by implementation of the modified Gibbs- 
Thomson relation at the interface: 


<tk 7V, (ci_ 

Tj ~ T '“ + P ,h i9 + I 


-Tc q ) 7 - 


(5) 


Jy 


pj{u, - uj) 


where 0 is the kinetic mobility. The above equation was first 
derived by Alexiades & Solomon (1993) and then modified 
by Juric & Tryggvason (1996) to include the kinetic mobility 
effect. 

The above equations are supplemented by the equations 
of state for the material properties: 


Vt 


where 
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= +u t ■ V. 
dt 


(7) 


Our numerical technique is an extension to the two-dimensi- 
onal model developed by Juric & Tryggvason (1995). Exten- 
sion of the method to three dimensions involves a number of 


complications in the operations on the front which will be ad- 
dressed in a future publication. We use a fixed (i.e. Eulerian), 
regular, and staggered grid to discretize the governing equa- 
tion and an unstructured (i.e. Lagrangian) grid to represent the 
front and to construct the material property fields. This un- 
structured grid is also used to distribute the source terms in the 
right hand side of mass conservation, momentum, and energy 
equations to the Eulerian grid and to interpolate the velocity, 
temperature, and density at the front points from the grid. 

The computations start with construction of the initial tem- 
perature field T n and the material property fields p , ft , c , k 
using known position of the front at / = 0, where n = 1 , The 
front is then moved to a new position using the interface veloc- 
ity at the current time TT t 71 . The surface tension is distributed 
to the grid as a body force using Peskin distribution function 
(see Peskin 1 977) and the density 1 and heat capacity r n+1 
fields are constructed at the new position. A modified projec- 
tion algorithm which is first order in time and second order in 
spatial dimensions is used to solve the momentum equation. 
In the absence of pressure term, the momentum equation is 
solved and a provisional (i.e. unprojected) velocity field u* is 
computed. The iterative part starts by guessing the front ve- 
locity at the next time step and construction of the mass source 
at the new position of the front. An elliptic pressure equa- 
tion is then solved and the provisional velocity is corrected 
for the pressure to obtain the projected velocity + 1 . The 
velocity of the fluid at the front position u j n+ 1 is interpolated 
using a Peskin interpolation function (sec Peskin 1977) and 
the heat source is constructed. The energy equation is then 
solved for T" +1 and the Gibbs-Thomson relation is checked 
at the front points. If this equation is satisfied for all the front 
points within a predetermined tolerance, the heat conduction 
coefficient k 7}+l and viscosity fields /i n+1 are updated for the 
next time step and the calculation proceeds. Otherwise, a new 
guess is proposed for the front velocity and the computations 
are repeated. 

The individual parameters that control the problem are 
pi , p v , pi, p V9 kiy k ( , O , c v , hjy, <7 , <!>, <t, T yo — T cq , and 
initial diameter of the bubble, do. Here, is the initial 
temperature of the fluid outside the bubble. Nondimcnsional- 
ization results the following parameters: p v j pi , Pv/pi , k v /ki, 
ejety Ja = ci(T.^ — T e(l )/hf y , Pr = pin/kiyGr = 
pig(pi — p v )do / pi 2 y C* o — uiTn Tt- q )(r / pih j g do, and 
i ) = ki/dpih jgdo. Ja is the Jacob number which is ratio 
of the sensible to the latent heat, Pr is the Prandtl number 
which is a measure of diffusion of momentum with respect to 
thermal energy, Gr has a strong resemblance to the Grashof 
number, Oa is a capillary number which is ratio of surface 
tension to viscous force, and 0 is nondimensional kinetic mo- 
bility. When we present our results, L = do, = pifpido, 
and t s = lju 6 are used as length, velocity, and time scales. 
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RESULTS 

Here, we present some preliminary results using our nu- 
merical technique. We are interested in growth of a bubble in 
a superheated liquid. The growth of a vapor bubble in a su- 
perheated liquid is controlled by three factors: liquid inertia, 
surface tension, and heat diffusion. These factors determine 
what is called “the critical radius” which is the minimum ra- 
dius that a nucleus should start with in order to grow. In the 
current study, we choose the initial bubble radius considerably 
larger than the critical radius. Therefore, we expect the bubble 
to grow. 

Wc start our analysis by considering the behavior of a sin- 
gle bubble in zero gravity. The first frame of figure (1) shows 
the initial position of the bubble in a domain. The domain 
size (nondimensionalized by dividing by bubble diameter) is 
2.5 x 2.5 x 5 and the grid resolution is 48 x 48 x 96. The 
initial front consists of 5832 triangular elements. Initially the 
liquid is superheated, the vapor is at the equilibrium tempera- 
ture, and the flow is quiescent. The nondimensional variables 
are la - 0.1, Pr = 0.25, Ca = 0.02, and 0 = 0.01. The 
ratio of material properties are p, .-//?/ = 0.1, p t /pi = 0.05, 
k L /k { = 0.025, and c v /c t = 1.0. The domain is periodic in 
the horizontal direction, wall -bounded (and insulated) at the 
bottom, and open at the lop. At the top boundary, the pressure 
and normal gradient of velocity and temperature are set to zero. 
The second frame shows the bubble and a normal section of 
velocity field at the middle of the domain at time / = 0.053. 
The bubble elongates in the vertical direction. Since the grav- 
ity is zero, there is no vortical structure inside the bubble. The 
velocity field in the vicinity of the bubble is relatively disturbed 
but it is uniform in the rest of the domain. The liquid exits from 
the top boundary due to fluid expansion at the interface. The 
fluid velocity is higher around the top of the bubble compared 
to its side. This is due to the constrain imposed on the flow 
by periodicity in the horizontal direction and the existence of 
the wall at the bottom. This results in a higher evaporation 
rate around the lop portion of the bubble and consequently the 
bubble deforms to an egg-like shape (third frame, / = 0.106). 
The ratio of final volume of the bubble to its initial volume is 
about fifteen. At the end of the simulation, the bubble consists 
of 45704 triangular elements. 

In the next simulation we study growth of a bubble un- 
der high gravity. The initial setup is the same as the previous 
simulation. The nondimensional numbers for this run are 
Ja = 0.1, Pr = 0.25, GV = 0.32 x 10\ Ca = 0.004, 
and 0 = 0.05. With the exception of density ratio which is 
increased to p v / pi = 0.5, all the other material property ratios 
arc the same as the corresponding ones in the first simulation. 
The grid resolution is the same as the previous run. Lower 
density difference results in a smaller front velocity which 
compensates for the higher resolution that may be needed to 


accurately resolve this case. The first frame of figure 2 shows 
the bubble and the velocity field at f = 0.0875. The velocity is 
higher inside the bubble due to buoyancy. Two counter rotat- 
ing vortices appear on the side of bubble which pump hot fluid 
from the top to the rear. Since the bubble grows, we expect the 
rise velocity to increase as a result of enhanced buoyancy. This 
is indeed the case as is seen from the scale of the velocities in 
the second frame at 1 = 0.22. The evaporation rate is more 
uniform compared to the previous case. This is due to the fact 
that the side of the bubble now has a chance to receive hot fluid 
as the bubble moves upward. As a result, the elongation in the 
vertical direction is much smaller compared to the first case. 
Moreover, unlike the previous run, an indentation appears at 
the rear of the bubble. This is a hydrodynamic effect and is 
formed as a result of local competition between the surface 
tension, pressure, and viscous forces. This indentation is in- 
creased as the bubble rises (third frame, i = 0.268). Although 
the high curvature developed at the rim of the bubble is a chal- 
lenge to accurate computation of surface tension, comparison 
of our results with curvatures of analytical surfaces showed 
that our method successfully resolves this issue. 

Inspection of the temperature field for both cases (not 
shown in the text) showed the gradual depletion of the super- 
heat. The initial sharp temperature gradient at the interface 
was replaced with a less steep temperature gradient as a result 
of heat diffusion. This did not have a pronounced effect at the 
early period of growth. However, it decreased the growth rate 
at the later time. 

CONCLUSIONS 

We have used direct numerical simulation to study phase 
change phenomena. Numerical simulation of a three-dimensio- 
nal boiling bubble in zero gravity showed the deformation of 
the bubble to an egg-like shape as a result of nonuniform evap- 
oration at the phase boundary. Similar compulations for a 
bubble under high gravity showed the formation of a “dim- 
pled” bubble. For both cases, the growth rate started to halt as 
a result of depletion of the superheat. 
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